Analiza danych krystalograficznych w R pochodzących z bazy Protein Data Bank (PDB)
Michał Hejduk
10 grudzień, 2018
Wstęp
Na początku pliku powinna znajdować się automatycznie wypełniona data generacji dokumentu oraz spis treści pozwalający przejść do najważniejszych sekcji. Ponadto raport powinien zaczynać się od rozdziału podsumowującego całą analizę, streszczającego najważniejsze spostrzeżenia analityka.
1. Wykorzystane biblioteki
library(dplyr)
library(tibble)
library(tidyr) #transformowanie danych
library(ggplot2)
library(plotly)
library(kableExtra) #ładne tabele
library(data.table) #%like%
library(corrplot)2. Powtarzalność wyników przy każdym uruchomieniu raportu na tych samych danych
set.seed(23)3.Wczytanie danych z pliku
Wczytujemy dane z pliku all_summary.csv. Za pomocą heurystyki wyznaczamy klasy pierwszych 100 wierszy.
filePath = "C:/Users/micha/Downloads/all_summary.csv"
initial <- read.table(filePath, nrows = 100, sep = ";", header = TRUE)
classes <- sapply(initial, class)
file_data <- read.table(filePath, colClasses = classes, sep = ";", na.strings = 'nan', nrows = 10000, header = TRUE)
data <- file_dataDane znajdują się w zmiennej data.
4.Wstępne czyszczenie danych:
Liczba wierszy w zbiorze: 10000.
Usuwamy wiersze posiadające wartość zmiennej ref_name równą: “UNK”, “UNX”, “UNL”, “DUM”, “N”, “BLOB”, “ALA”, “ARG”, “ASN”, “ASP”, “CYS”, “GLN”, “GLU”, “GLY”, “HIS”, “ILE”, “LEU”, “LYS”, “MET”, “MSE”, “PHE”, “PRO”, “SEC”, “SER”, “THR”, “TRP”, “TYR”, “VAL”, “DA”, “DG”, “DT”, “DC”, “DU”, “A”, “G”, “T”, “C”, “U”, “HOH”, “H20”, “WAT”
res_names_to_remove <- c('UNK', 'UNX', 'UNL', 'DUM', 'N', 'BLOB', 'ALA', 'ARG', 'ASN', 'ASP', 'CYS', 'GLN', 'GLU', 'GLY', 'HIS', 'ILE', 'LEU', 'LYS', 'MET', 'MSE', 'PHE', 'PRO', 'SEC', 'SER', 'THR', 'TRP', 'TYR', 'VAL', 'DA', 'DG', 'DT', 'DC', 'DU', 'A', 'G', 'T', 'C', 'U', 'HOH', 'H20', 'WAT')
data <- filter(data, !res_name %in% res_names_to_remove)Liczba wierszy w zbiorze po usunięciu wierszy z ww. wartościami res_name: 9940.
5. Analiza brakujących danych:
W poniższej sekcji oczyścimy dane z kolumn i wierszy zawierających wartości NaN.
Na początku analiza kolumn.
Wyznaczamy kolumny, które zawierają przynajmniej jedną wartość NaN, grupujemy po nazwie.
column_value_DT <- gather(data, "column_name", "val", 1:ncol(data))
columns_with_nans <- filter(column_value_DT, is.na(column_value_DT$val))[,1]Wyświetlamy podsumowanie:
table(columns_with_nans)## columns_with_nans
## dict_atom_C_count dict_atom_N_count
## 174 174
## dict_atom_non_h_count dict_atom_non_h_electron_sum
## 174 174
## dict_atom_O_count dict_atom_S_count
## 174 174
## part_01_density_CI part_01_density_E2_E1
## 80 80
## part_01_density_E3_E1 part_01_density_E3_E2
## 80 80
## part_01_density_FL part_01_density_FL_norm
## 80 80
## part_01_density_I1 part_01_density_I1_norm
## 80 80
## part_01_density_I2 part_01_density_I2_norm
## 80 80
## part_01_density_I3 part_01_density_I3_norm
## 80 80
## part_01_density_I4 part_01_density_I4_norm
## 80 80
## part_01_density_I5 part_01_density_I5_norm
## 80 80
## part_01_density_I6 part_01_density_I6_norm
## 80 80
## part_01_density_M000 part_01_density_O3
## 80 80
## part_01_density_O3_norm part_01_density_O4
## 80 80
## part_01_density_O4_norm part_01_density_O5
## 80 80
## part_01_density_O5_norm part_01_density_sqrt_E1
## 80 80
## part_01_density_sqrt_E2 part_01_density_sqrt_E3
## 80 80
## part_01_density_Z_0_0 part_01_density_Z_1_0
## 80 80
## part_01_density_Z_2_0 part_01_density_Z_2_1
## 80 80
## part_01_density_Z_3_0 part_01_density_Z_3_1
## 80 80
## part_01_density_Z_4_0 part_01_density_Z_4_1
## 80 80
## part_01_density_Z_4_2 part_01_density_Z_5_0
## 80 80
## part_01_density_Z_5_1 part_01_density_Z_5_2
## 80 80
## part_01_density_Z_6_0 part_01_density_Z_6_1
## 80 80
## part_01_density_Z_6_2 part_01_density_Z_6_3
## 80 80
## part_01_density_Z_7_0 part_01_density_Z_7_1
## 80 80
## part_01_density_Z_7_2 part_01_density_Z_7_3
## 80 80
## part_01_shape_CI part_01_shape_E2_E1
## 80 80
## part_01_shape_E3_E1 part_01_shape_E3_E2
## 80 80
## part_01_shape_FL part_01_shape_FL_norm
## 80 80
## part_01_shape_I1 part_01_shape_I1_norm
## 80 80
## part_01_shape_I2 part_01_shape_I2_norm
## 80 80
## part_01_shape_I3 part_01_shape_I3_norm
## 80 80
## part_01_shape_I4 part_01_shape_I4_norm
## 80 80
## part_01_shape_I5 part_01_shape_I5_norm
## 80 80
## part_01_shape_I6 part_01_shape_I6_norm
## 80 80
## part_01_shape_M000 part_01_shape_O3
## 80 80
## part_01_shape_O3_norm part_01_shape_O4
## 80 80
## part_01_shape_O4_norm part_01_shape_O5
## 80 80
## part_01_shape_O5_norm part_01_shape_sqrt_E1
## 80 80
## part_01_shape_sqrt_E2 part_01_shape_sqrt_E3
## 80 80
## part_01_shape_Z_0_0 part_01_shape_Z_1_0
## 80 80
## part_01_shape_Z_2_0 part_01_shape_Z_2_1
## 80 80
## part_01_shape_Z_3_0 part_01_shape_Z_3_1
## 80 80
## part_01_shape_Z_4_0 part_01_shape_Z_4_1
## 80 80
## part_01_shape_Z_4_2 part_01_shape_Z_5_0
## 80 80
## part_01_shape_Z_5_1 part_01_shape_Z_5_2
## 80 80
## part_01_shape_Z_6_0 part_01_shape_Z_6_1
## 80 80
## part_01_shape_Z_6_2 part_01_shape_Z_6_3
## 80 80
## part_01_shape_Z_7_0 part_01_shape_Z_7_1
## 80 80
## part_01_shape_Z_7_2 part_01_shape_Z_7_3
## 80 80
## part_02_density_CI part_02_density_E2_E1
## 704 704
## part_02_density_E3_E1 part_02_density_E3_E2
## 704 704
## part_02_density_FL part_02_density_FL_norm
## 704 704
## part_02_density_I1 part_02_density_I1_norm
## 704 704
## part_02_density_I2 part_02_density_I2_norm
## 704 704
## part_02_density_I3 part_02_density_I3_norm
## 704 704
## part_02_density_I4 part_02_density_I4_norm
## 704 704
## part_02_density_I5 part_02_density_I5_norm
## 704 704
## part_02_density_I6 part_02_density_I6_norm
## 704 704
## part_02_density_M000 part_02_density_O3
## 704 704
## part_02_density_O3_norm part_02_density_O4
## 704 704
## part_02_density_O4_norm part_02_density_O5
## 704 704
## part_02_density_O5_norm part_02_density_sqrt_E1
## 704 704
## part_02_density_sqrt_E2 part_02_density_sqrt_E3
## 704 704
## part_02_density_Z_0_0 part_02_density_Z_1_0
## 704 704
## part_02_density_Z_2_0 part_02_density_Z_2_1
## 704 704
## part_02_density_Z_3_0 part_02_density_Z_3_1
## 704 704
## part_02_density_Z_4_0 part_02_density_Z_4_1
## 704 704
## part_02_density_Z_4_2 part_02_density_Z_5_0
## 704 704
## part_02_density_Z_5_1 part_02_density_Z_5_2
## 704 704
## part_02_density_Z_6_0 part_02_density_Z_6_1
## 704 704
## part_02_density_Z_6_2 part_02_density_Z_6_3
## 704 704
## part_02_density_Z_7_0 part_02_density_Z_7_1
## 704 704
## part_02_density_Z_7_2 part_02_density_Z_7_3
## 704 704
## part_02_shape_CI part_02_shape_E2_E1
## 704 704
## part_02_shape_E3_E1 part_02_shape_E3_E2
## 704 704
## part_02_shape_FL part_02_shape_FL_norm
## 704 704
## part_02_shape_I1 part_02_shape_I1_norm
## 704 704
## part_02_shape_I2 part_02_shape_I2_norm
## 704 704
## part_02_shape_I3 part_02_shape_I3_norm
## 704 704
## part_02_shape_I4 part_02_shape_I4_norm
## 704 704
## part_02_shape_I5 part_02_shape_I5_norm
## 704 704
## part_02_shape_I6 part_02_shape_I6_norm
## 704 704
## part_02_shape_M000 part_02_shape_O3
## 704 704
## part_02_shape_O3_norm part_02_shape_O4
## 704 704
## part_02_shape_O4_norm part_02_shape_O5
## 704 704
## part_02_shape_O5_norm part_02_shape_sqrt_E1
## 704 704
## part_02_shape_sqrt_E2 part_02_shape_sqrt_E3
## 704 704
## part_02_shape_Z_0_0 part_02_shape_Z_1_0
## 704 704
## part_02_shape_Z_2_0 part_02_shape_Z_2_1
## 704 704
## part_02_shape_Z_3_0 part_02_shape_Z_3_1
## 704 704
## part_02_shape_Z_4_0 part_02_shape_Z_4_1
## 704 704
## part_02_shape_Z_4_2 part_02_shape_Z_5_0
## 704 704
## part_02_shape_Z_5_1 part_02_shape_Z_5_2
## 704 704
## part_02_shape_Z_6_0 part_02_shape_Z_6_1
## 704 704
## part_02_shape_Z_6_2 part_02_shape_Z_6_3
## 704 704
## part_02_shape_Z_7_0 part_02_shape_Z_7_1
## 704 704
## part_02_shape_Z_7_2 part_02_shape_Z_7_3
## 704 704
## weight_col
## 9940Jak widać na podsumowaniu powyżej, dla wszystkich wierszy kolumna weight_col ma wartość NaN. Usuwamy ją ze zbioru danych:
data <- data[,-which(names(data) == 'weight_col' )]Widać również, że część wierszy nie posiada wartości dla kolumn zaczynających się od part_02 oraz part_01.
W celu oczyszczenia danych usuwamy te wiersze.
all_rows_count <- nrow(data)
data <- data[complete.cases(data),]
cleaned_rows_count <- nrow(data)
removed_rows_count <- all_rows_count - cleaned_rows_count Usunięto 869 wierszy. Stanowi to 8.7424547 % wszystkich wierszy.
6. Podsumowanie danych:
Liczba wierszy: 9071
Liczba kolumn: 411
Liczba unikalnych wartości res_name: 801
Poniżej rozmiar data frame w pamięci (w Mb):
print(object.size(data), units='Mb')## 36.8 Mb7. Ograniczenie liczby klas do 50 najpopularniejszych wartości
W tej sekcji filtrujemy dane do wierszy, których klasa znajduje się wśród 50 najpopularniejszych.
Na początku grupujemy dane po atrybucie res_name i sortujemy.
top50 <- sort(table(data$res_name), decreasing=TRUE)[1:50]
top50names <- names(top50)Poniżej znajduje się posortowana lista 50 najpopularniejszych klas.
top50names## [1] "SO4" "GOL" "EDO" "NAG" "CL" "DMS" "ZN" "CA" "HEM" "MG" "PO4"
## [12] "NA" "NAD" "ACT" "FAD" "PEG" "IOD" "CD" "K" "PG4" "MN" "NAP"
## [23] "FMN" "MLY" "LHG" "AMP" "NDP" "FE2" "EPE" "ADP" "PLP" "ACY" "HEC"
## [34] "BR" "MAN" "CU" "FE" "FMT" "MPD" "C8E" "CIT" "CME" "PGE" "FES"
## [45] "H4B" "PLC" "CYC" "FEC" "GLC" "NI"Filtrujemy za pomocą wyznaczonych 50 nazw
data <- filter(data, res_name %in% top50names)
data$res_name <- droplevels(data$res_name)Pozostało 6376 wierszy.
8.Korelacja między zmiennymi:
Korelację wyznaczamy dla kolumn numerycznych.
Pomijamy kolumny: “blob_coverage” “res_coverage” “title” “pdb_code” “res_name” “chain_id” “skeleton_data” “fo_col” “fc_col”
Dla zbioru danych korelację wyznaczamy za pomocą funkcji cor.
Następnie prezentujemy wyniki, za pomocą funkcji corrplot.
is_numeric_column_list <- unlist(lapply(data, is.numeric))
numeric_data <- data[, is_numeric_column_list]
cor_data <- cor(numeric_data, use = "pairwise.complete.obs")
#pdf(file = "korelacja_miedzy_kolumnami_FULL.pdf")
corrplot(cor_data, method = "color", tl.col="black",
tl.cex = 0.2, cl.cex = 0.2, title = "Korelacja między kolumnami")
#dev.off()Z racji dużego rozmiaru macierzy (402 x 402) wykres jest nieczytelny.
Wynik można zobaczyć również w wyeksportowanym pliku pdf: korelacja_miedzy_kolumnami.pdf
W celu zwiększenia widoczności, ogarniczymy liczbę kolumn, tylko do tych zaczynających się od part_01:
is_part_01_column <- colnames(numeric_data) %like% 'part_01'
part_01_column_names <- colnames(numeric_data)[is_part_01_column]
filtered_numeric_data <- numeric_data[is_part_01_column]
cor_data2 <- cor(filtered_numeric_data, use = "pairwise.complete.obs")
#pdf(file = "korelacja_miedzy_kolumnami_part_01.pdf")
corrplot(cor_data2, tl.col="black", order = "FPC",
tl.cex = 0.3, title = "Korelacja między kolumnami")
#dev.off()W pliku “korelacja_miedzy_kolumnami_part_01.pdf” znajduje się powyższa wizualizacja.
9. Liczba przykładów dla każdej klasy
sort(table(data$res_name), decreasing=TRUE)##
## SO4 GOL EDO NAG CL DMS ZN CA HEM MG PO4 NA NAD ACT FAD PEG IOD CD
## 982 587 466 396 380 327 318 282 234 194 169 161 129 116 112 104 96 87
## K PG4 MN NAP FMN MLY LHG AMP NDP FE2 EPE ADP PLP ACY HEC BR MAN CU
## 76 71 59 55 54 54 52 45 45 43 42 41 41 37 36 33 33 32
## FE FMT MPD C8E CIT CME PGE FES H4B PLC CYC FEC GLC NI
## 32 31 29 28 28 28 28 27 27 27 26 26 25 2510. Wykresy rozkładów liczby atomów (local_res_atom_non_h_count) i elektronów (local_res_atom_non_h_electron_sum).
chart_data_1 <- ggplot(data, aes(x = data$local_res_atom_non_h_count))
chart_layout_1 <- geom_histogram(col = 'blue', fill = 'orange', binwidth = 1)
chart1 <- chart_data_1 + chart_layout_1 + labs(title="Rozkład liczby atomów") + labs(x="Liczba atomów")
chart1
chart_data_2 <- ggplot(data, aes(x = data$local_res_atom_non_h_electron_sum))
#chart_layout_2 <- geom_histogram(col = 'blue', fill = 'orange', binwidth = 5)
chart2 <- chart_data_2 + chart_layout_1 + labs(title="Rozkład liczby elektonów") + labs(x="Liczba elektronów")
chart2
11. Analiza klas pod względem niezgodności liczby atomów i elektronów
Na początek wybieramy interesujące nas kolumny do porównań: local_res_atom_non_h_count, oraz dict_atom_non_h_count. Dodatkowo każdy wiersz zawiera nazwę klasy.
atoms_incompatibility <- unique(select(data, res_name, local_res_atom_non_h_count, dict_atom_non_h_count))Dodajemy kolumnę z różnicą, wyliczoną na postawie 2 innych.
atoms_incompatibility_with_diff <- mutate(atoms_incompatibility, diff = abs(local_res_atom_non_h_count - dict_atom_non_h_count))Sortujemy i wyświetlamy 5 pierwszych.
atoms_incompatibility_with_diff_sorted <- atoms_incompatibility_with_diff[order(-atoms_incompatibility_with_diff$diff),]
head(atoms_incompatibility_with_diff_sorted, 5)## res_name local_res_atom_non_h_count dict_atom_non_h_count diff
## 63 LHG 8 49 41
## 62 LHG 11 49 38
## 77 PLC 7 42 35
## 78 PLC 9 42 33
## 4 NDP 22 48 26Poniżej znajduje się tabela pokazująca 10 klas z największą niezgodnością liczby atomów:
atoms_table_data <- select(atoms_incompatibility_with_diff_sorted, res_name, diff)[1:10,]
row.names(atoms_table_data) <- NULL
atoms_table_data %>% knitr::kable(col.names = c('res_name', 'Niezgodność liczy atomów')) %>% kable_styling()| res_name | Niezgodność liczy atomów |
|---|---|
| LHG | 41 |
| LHG | 38 |
| PLC | 35 |
| PLC | 33 |
| NDP | 26 |
| PLC | 24 |
| PLC | 23 |
| NAP | 22 |
| PLC | 22 |
| NAP | 17 |
Analogicznie postępujemy aby pokazać niezgodność liczby elektronów:
electrons_incompatibility <- unique(select(data, res_name, local_res_atom_non_h_electron_sum, dict_atom_non_h_electron_sum))
electrons_incompatibility_with_diff <- mutate(electrons_incompatibility, diff = abs(local_res_atom_non_h_electron_sum - dict_atom_non_h_electron_sum))
electrons_incompatibility_with_diff_sorted <- electrons_incompatibility_with_diff[order(-electrons_incompatibility_with_diff$diff),]
electrons_table_data <- select(electrons_incompatibility_with_diff_sorted, res_name, diff)[1:10,]
row.names(electrons_table_data) <- NULL
electrons_table_data %>% knitr::kable(col.names = c('res_name', 'Niezgodność liczy elektronów')) %>% kable_styling()| res_name | Niezgodność liczy elektronów |
|---|---|
| LHG | 275 |
| LHG | 257 |
| PLC | 236 |
| PLC | 224 |
| NDP | 170 |
| PLC | 170 |
| PLC | 164 |
| NAP | 158 |
| PLC | 158 |
| NAP | 112 |
12. Rozkład wartości wszystkich kolumn zaczynających się od part_01
Na początek wyszukujemy nazwy kolumn zaczynające się od part_01. Wyświetlamy 10 pierwszych.
is_part_01_column <- colnames(data) %like% 'part_01'
part_01_column_names <- colnames(data)[is_part_01_column]
head(part_01_column_names, 10)## [1] "part_01_shape_segments_count" "part_01_density_segments_count"
## [3] "part_01_volume" "part_01_electrons"
## [5] "part_01_mean" "part_01_std"
## [7] "part_01_max" "part_01_max_over_std"
## [9] "part_01_skewness" "part_01_parts"Następnie filtrujemy dane, przekształcamy do postaci:
filtered_data <- data[is_part_01_column]
gathered_data <- gather(filtered_data, 'col_name', "value", part_01_column_names)
dim(gathered_data)## [1] 675856 2wyświetlamy rozkład wartości, osobny wykres dla każdej kolumny:
k <- ggplot(gathered_data, aes(x = gathered_data$value))
i <- k + facet_wrap(vars(gathered_data$col_name), ncol = 3, scales = "free")
i + geom_histogram(bins = 25) + xlab("Rozkład wartości")
13.Interaktywny wykres z liczbą atomów, elektronów, oraz nazwą klasy:
p <- ggplot(data, aes(local_res_atom_non_h_count, local_res_atom_non_h_electron_sum,
color=res_name)) +
geom_point()
ggplotly(p)5.Brakujące dane:
5.Brakujące dane:
Including Plots
You can also embed plots, for example:
Note that the echo = FALSE parameter was added to the code chunk to prevent printing of the R code that generated the plot.